Fluctuation-Dissipation Relations for Motions of Center of Mass in Driven Granular 

Fluids under Gravity 
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We investigate a validity of fluctuation-dissipation relations in a nonequilibrium stationary state 
of fluidized granular media under gravity. A phenomenological Langevin-type theory describing the 
fluctuation of center of mass height, which was originally constructed for one-dimensional granular 
gas on a vibrating bottom plate, is generalized to any dimensions even for the case that the vibrating 
bottom plate is replaced by a thermal wall. The theory gives analytical expressions for the power 
spectrum and response function of the center of mass height; furthermore, it predicts a fluctuation- 
dissipation relation between them, which is known to be satisfied at equilibrium, with a modification 
that equilibrium temperature is replaced by an effective temperature defined by a kinetic energy 
of the center of mass. To check these explicit theoretical predictions, we performed extensive and 
accurate event-driven molecular dynamics simulations for the model system with a thermal wall at 
the bottom. We found that the power spectrum and response function of the center of mass height 
show good agreement with theoretical predictions within a range of time scales in which our theory 
is valid. As the most remarkable result, it is shown that a fluctuation-dissipation relation for the 
granular system is well satisfied especially at a large frequency (short time) region in a wide range 
of system parameter. We finally remark that the relation between systematic deviations at a small 
frequency (long time) region and time scales of the driven granular system. 

PACS numbers: 45.70.-n, 47.70.Nd, 05.40.-a, 05.70.Ln 



I. INTRODUCTION 

Granular materials show fluid like behavior when they 
are supplied sufficient energy by external vibration. Such 
fluidized states of granular matter have been studied as 
an interesting example of nonequilibrium fluids that ex- 
hibits a rich variety of phenomena such as convection, 
pattern formation on the surface, and segregation (see 
Ref. [l| and references therein). Besides these pattern 
forming instabilities, a plain stationary state of vibrated 
granular fluids without any complex spatial structures 
serves as an archetypal example of nonequilibrium sta- 
tionary states (NESSs). It has been a fundamental sub- 
ject for long years to find any thcrmodynamic-like de- 
scription or to identify common property of fluctuations 
in a wide variety of NESS in nature. 

One of important issues, addressed in this paper, is 
to clarify the validity of a fluctuation-dissipation relation 
(FDR) in granular fluids subject to external vibration. A 
FDR connects the response of an equilibrium system to 
a small perturbation with the time correlation of spon- 
taneous fluctuations in the system without perturbation. 
It has attracted recently a lot of interest how a FDR is 
violated or should be modified in ageing systems such as 
glass and in NESSs of various systems. (See Refs. [H, H| 
as recent reviews.) 

For granular systems, FDRs have been studied for sev- 
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eral situations. Many works have devoted to the case of 
freely cooling granular gas, in which the gas develops 
freely without external forces and "cools" down as a re- 
sult of dissipative nature of the grain interactions, aiming 
to derive a (modified) Grccn-Kubo relation and to calcu- 
late transport coefficients using it 0-0|- While there is 
no stationary state of freely cooling granular g NESS 
of granular gas can be achieved by supplying energy from 
outside by means of external driving. A typical exper- 
imental mean of injecting energy is shaking a container 
or vibrating a bottom wall (sec, e.g., Ref. Q). In the 
case that the shaking (or vibrating) is strong enough to 
inject energy to all grains by frequent collisions with the 
vibrating wall, the effect of the vibrating wall is often 
modeled using a thermal bath that couples to every par- 
ticle. FDRs in such uniforml y d riven granular systems 
have been studied in Refs. [9l-tl3|. Puglisi et al. Q car- 
ried out numerical simulations of a model of uniformly 
driven granular gas and studied FDRs for two different 
observables. They observed the FDRs are well satisfied 
if the equilibrium temperature in the FDR for a system 
at equilibrium is replaced by the granular temperature 
defined as the mean-square fluctuation of the grain ve- 
locity. Garzo [Io[ studied analytically diffusion of im- 
purities immersed in a granular gas under influence of 
uniform driving forces and showed that a modified form 
of the Einstein relation, in which the temperature of the 
gas is replaced by the temperature of the impurity, is vio- 
lated due to the non-Maxwellian behavior of the impurity 
velocity distribution function. Bunin et al. [l2j analyzed 
a mean-field model of uniformly driven granular gas and 
showed that an effective temperature defined by a FDR 
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depends on frequency. In the case that the shaking (or 
vibrating) is not strong enough to be regarded as uni- 
form driving, energy injection through a boundary has 
to be explicitly considered. Brey et al. [l3| studied nu- 
merically the fluctuations of volume of a vibrated low 
density granular gas confined by a mobile piston on the 
top. In this system energy is supplied from the vibrating 
bottom wall. They discussed interpretation of an effec- 
tive temperature defined by requiring the same relation 
between fluctuations of volume and compressibility as in 
equilibrium systems. FDRs and effective temperatures 
in much denser systems have also been studied by sev- 
eral authors [lj-[l7|. Among these studies, we refer to 
an experimental study by D'Anna et al. |17| . since their 
theory based on a Langevin equation formally has the 
same form as ours although their experimental setup is 
very different from ours. They carried out an experiment 
observing fluctuating motion of a torsion oscillator im- 
mersed in vibration-fluidized granular matter and found 
that it can be described as a first approximation by a 
formalism on Brownian motion in the equilibrium and a 
FDR with an effective temperature approximately holds. 

This paper is intended to investigate fluctuating mo- 
tion of the center of mass (COM) in a NESS of a granular 
matter fluidized by an external energy source located at a 
bottom wall, under influence of gravity. Instead of using 
macroscopic probes such as a piston fl3j and a torsion 
oscillator [l?} , we focused on position of the COM that is 
observable by means of digital high speed photography in 
experiments (l8j . Our major motivation to study fluctua- 
tions of the COM is that a simple (or universal) law might 
hold as a result of the following properties: First, the 
fluctuations of macrovariables like position of the COM 
often possess the largest time scale in the system. Sec- 
ondly, they are expected to have a Gaussian property in 
a similar sense as the central limit theorem. (In the case 
of a Markovian stochastic process, a Gaussian property 
of fluctuations of macrovariables can indeed be derived 
from a master equation of the Markovian process [19|.) 
With this expectation, we proposed a phcnomenological 
theory based on a simple formalism on Brownian motion 
that describes motion of the COM height in a NESS of 
a one-dimensional vibrated granular fluid (20j ; we found 
that the important qualitative features of dynamics of 
the COM in event-driven molecular dynamics (MD) sim- 
ulations are all accounted for by the theory. The theory 
was extended to a two-dimensional granular fluid on a 
thermal wall [2l|. In this paper we will show that when 
we apply the phcnomenological theory to granular flu- 
ids in higher dimensions, careful consideration of time 
scales in granular hydrodynamics (22l . |23| is necessary. 
Within a time range in which our theory is valid, it pre- 
dicts an existence of a FDR except that the equilibrium 
temperature in the FDR for an equilibrium system must 
be modified by an effective temperature of velocity fluc- 
tuation of the COM. In order to test our prediction, we 
performed extensive and accurate event-driven MD sim- 
ulations for a two-dimensional system of inelastic hard 



disks on a thermal wall. 

Our main result is that a FDR with an effective tem- 
perature holds within statistical uncertainty in simula- 
tions in a large frequency (short time) region while it is 
violated in a small frequency (long time) region. The 
effective temperature is defined by a kinetic energy of 
the COM. We observed in our simulations that the ratio 
between the effective temperature and the global gran- 
ular temperature increases with inelasticity; the former 
can be more than 4 times larger than the latter for the 
highest inelasticity case. 

This paper is organized as follows. In Sec. II, we de- 
scribe a model granular system and discuss important 
time scales in the system. In Sec. Ill, the Langevin equa- 
tion is introduced, and analytical expressions for power 
spectrum and response function of the COM height are 
expressed briefly. We also remark on the FDR between 
these two functions. Note that a complete derivation 
of the Langevin equation and detail calculation for the 
power spectrum and the response function are summa- 
rized in Appendix A and B, respectively. As the numeri- 
cal tests, the comparison between the theoretical predic- 
tions and extensive event-driven MD are shown in Sec. IV. 
Finally, in Sec.V, we summarize the main results of a va- 
lidity of FDR in this paper and remark on the relation 
between systematic deviations at a small frequency (long 
time) region and time scales of the driven granular sys- 
tem. 



II. THE MODEL SYSTEM 
A. System 

As a model of grains bouncing on a vibrating bottom 
plate under gravity, we consider a d-dimensional system 
of N inelastic particles on a "thermal" bottom wall in a 
constant gravitational field g. The particles in the system 
have diameter a and mass m; the total mass of particles 
is denoted by M (= Nm). The thermal wall is kept 
at a constant temperature To, which plays a role of a 
heat source supplying sufficient translational energy to 
the particles to fluidize them. The z-direction is chosen 
to be opposite to the direction of gravity, and the thermal 
wall is fixed at z = 0. For simplicity, we adopt periodic 
boundary conditions in horizontal directions in order to 
ignore the boundary (side-wall) effects. Collisions be- 
tween particles are inelastic; inelasticity of the particle 
collisions is characterized by a normal restitution coeffi- 
cient r. In order to avoid any pattern forming instability 
along horizontal directions, we choose both the inelas- 
ticity and linear scales of the system in horizontal direc- 
tions are sufficiently small so that the system remains 
homogeneous in horizontal directions. These conditions 
of setting will be discussed in more detail in Sec. IV. 
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B. Time scales 

Before discussing on important time scales in the 
system, let us define several quantities that charac- 
terize macroscopic property of the system. We first 
define the kinetic energy per particle as K(t) = 
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i(t) 2 /2, the long time average of Kit) in 



a NESS as K = lim T _, 00 (l/T) / Q K(t)dt (hereafter the 
ovcrline on a quantity represents the long time average of 
the quantity in a NESS), the global granular temperature 
T by ksT = (2/d)K, where fc^ is the Boltzmann con- 
stant, and the thermal velocity as c = (rife^T '/to) 1 / 2 = 
(2K /m) 1 / 2 . A characteristic length scale of the system 
in vertical direction I is then defined as / = c 2 /g. 

Bromberg et al. [23| have suggested that there are three 
important time scales in this system at a level of hydro- 
dynamics: macroscopic oscillation time t osc (referred to 
as "fast time scale" in Ref . [23| ) , relaxation time for ther- 
mal conduction Ttherm > an d relaxation time for collisional 
dissipation r<jj ss . 

For the sake of simplicity in estimation of these time 
scales, let us assume here the system is nearly homoge- 
neous although this is not true for small r and large N. 
The time scale t osc represents a period of the slowest os- 
cillation in vertical direction, that is, a period of sound 
mode with the longest wavelength. Thus t osc ~ l/c s , 
where c s is a sound velocity. Assuming c a ~ c which is 
well satisfied in the case of normal gas, r osc can be es- 
timated as t osc r~j c/g. Since l/c s also characterize the 
relaxation time of pressure, t p , we can regard t p and t osc 
are on the same order: r p ~ t osc ~ c/g. The relax- 
ation time for thermal conduction Ttherm, is estimated as 
Ttherm ~ ^ 2 /( k //° c p)j where k is the thermal conductivity, 
p is the mass density, and c p is the specific heat at con- 
stant pressure j24[. p can be estimated as p ~ Mj (I A) ~ 
mN z /{la d ~ 1 ), where A represents the area of the bot- 
tom plate in the case of three dimensions (A represents 
the length of the bottom plate in two dimensions and 
A = 1 in one dimension) and N z represents the number 
of monolayers at rest, k and c p are obtained from kinetic 
theory for elastic spheres and disks [25j]: k ~ kBc/a d ~ 1 
and c p ~ fcs/m. Substituting these results, we obtain 
Ttherm ~ N z c/g. The relaxation time for collisional dis- 
sipation Tdiss can be estimated as the inverse of (1 — r 2 )v, 
where v is the collision frequency between two parti- 
cles. Substituting the lowest order estimation of v based 
on kinetic theory, v ~ pcr d ~ 1 c/m ~ N z c/l, we obtain 
T dlss ~ (1 - r^c/iNzg). 

The time scales estimated above are summarized as 
follows: 



that for a system with given N z and r, macroscopic dy- 
namics with time that is scaled by c/g are independent 
of g. We will utilize this fact later in order to obtain a 
frequency response function in an efficient way. 

There are three dimcnsionless parameters obtained 
from these three time scales. The first is T t herm/ T osc ~ 
N z . The second is t osc /t,hss ~ N z (l — r 2 ). The third is 
Ttherm / Tdiss ~ N 2 {1 — r 2 ) . The first and third parame- 
ters are governing parameters of hydrodynamic descrip- 
tion of the system introduced by Bromberg et al. [231 ] . 
They have shown that steady state profile is governed 
only by a parameter 



A^iV^l-r 2 ) 1 / 2 , 



(2) 



which is proportional to (Ttherm / 1 Tdiss) 1 ^ 2 ■ If r <C 1, the 
second parameter T osc /Tdi SS is related to X = N z (l — r) 
that plays a role of a governing parameter of transition 
from a condensed phase into a fluidized phase in a system 
of one-dimensional column of beads on vibrating bottom 
plate |26j . In our study, we consider the case N z ^> 
land assume Ttherm S> T osc i n the following theoretical 
analysis. 



III. THEORETICAL DERIVATIONS OF 
FLUCTUATION-DISSIPATION RELATION 

In this section, we summarize the theoretical deriva- 
tion of (i) the power spectrum, (ii) frequency response 
function, and (iii) FDR between (i) and (ii). Firstly, we 
introduce a Langevin equation as a first approximation 
which describes fluctuating motion of the COM within 
the fast time scale of r osc and t p . Note that the deriva- 
tion of our theory has already been published in Ref. [2(| ■ 
We assume Ttherm 3> t osc as mentioned above, and focus 
on dynamics of the COM in the time scale of t osc ignor- 
ing a significant slow relaxation process of fluctuations of 
global granular temperature around its stationary value 
(2/ d)K /ks- The effect of this slow dynamics of granular 
temperature and validity of our assumption of time scale 
will be discussed later. 

We summarize details of derivation of our Langevin 
formalism in Appendix A and show here the final result. 
Let us denote the height of the COM of granular fluids 
at time t as Z(t), a time average of Z(t) over a long time 
interval in a NESS as Z, and small deviation of Z(t) 
from Z as SZ(t) = Z(t) — Z. A Langevin equation for 
fluctuating motion of 8Z(t) is given by (see Eq. (|A5I) in 
Appendix A), 
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It is important to note that the all time scales, r osc , t„, 
Ttherm, and Tdiss oxe proportional to c/g. This means 



d 2 5Z 2 d6Z 
, „ = -il oZ — u—r~ 
dt 2 ^ dt 



R(t) 

M ' 



(3) 



where R(t) represents a random force which is assumed 
to be a Gaussian white noise: 



(R(t)) = 0, (R{t)R(t r )) = IS(t - t'). 



(4) 
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The brackets 



denote an average over the random 



si 



force. In NESS, it is reasonable to assume (Z(t)) 
where (• • • ) st represents the average in a stationary state. 
The constant / represents intensity of the random force, 
which is related to the second moment of the velocity 
fluctuations of the COM. This relation can be obtained 
by calculating the average kinetic energy of the COM mo- 
tion in z-direction per particle: Kqmz = (MV z (t) 2 /2) , 
here V z is the z-component of the velocity of the COM: 
V z (t) = dZ }f^ ■ Using an analytical solution Eq. (|B1[) 
of the Langevin equation, we obtain Kcmz = I/AMfi. 
Hence, the constant I is identified as follows: 



4M>if, 



CMz- 



(5) 



This is indeed the same procedure to determine the 
noise intensity / when the Langevin equation describes 
fluctuations in equilibrium at temperature T. In equi- 
librium, equipartition of energy implies Kqmz — ksT/2. 
Thus, we obtain the well-known result I = 2M [iksT. In 
the case of NESS of granular fluids, violation of equiparti- 
tion of energy is observed in various systems (see Ref. [27j 
and references therein). We will show later a result 
of numerical simulations that clearly shows violation of 
equipartition, Kcmz 7^ k B T /2, when we recognize T as 
the global granular temperature. 

The coefficients f2 and fi describe an angular frequency 
of the slowest oscillation of the COM height and frictional 
coefficient with respect to relative motion of the COM 
height against the bottom wall, respectively. According 
to the time scales we consider here, we assume Q ~ r, 



and /i 



and write them as, 



Q = fig I c, fx = jig/ c. 
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(6) 



Since values of the coefficients Q and jl cannot be esti- 
mated in our phenomcnological theory, they are fixed as 
fitting parameters when we compare results of simula- 
tions with the theoretical predictions. 

Power Spectrum: The power spectrum S(oj) that 
represents property of fluctuations of Z around NESS is 
defined as the Fourier transform of the time correlation 
function, 



dt e 



(SZ(0)SZ(t)) s 



(7) 



Derivation of S(u>) using the analytic solution of the 
Langevin equation is straightforward. We express here 
the final expressions of S(ui) in this system (Sec Appendix 
B for more detail deviation), 



S(u) = 



1 



CMz 



M (ft2 
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Response Function: The frequency response func- 
tion which characterizes linear response of Z in 
NESS against a small external force ef(t) can be defined 



as, 



x(uj) = lim 

e— ¥0 



SZ(oj) )/ef(u), 



(9) 



where SZ(u>) and f(u>) are the Fourier transform of SZ(t) 
and /(i), respectively. The analytical expression of x(w) 
is given by (See Appendix B for derivation), 
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(10) 



According to conventional definition, x( w ) can be de- 
composed into real x'( w ) an d imaginary x"( w ) parts as 
x(w) = x'i^) ~ *x"( CJ )- Thus we obtain the following 
expressions, 



x'M 



f 
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M(ft 2 -w 2 ) 2 + (Aiw) 2 ' 



X» = 
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~M (fl 2 — uj 2 ) 2 + ([ioj)' 



(11) 
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Fluctuation Dissipation Relation: Comparing 
Eqs. © and ([T2")l. we obtain a following FDR, 



2k B T eff 



X», 



(13) 



where T e f f is an effective temperature defined as T e f f = 
2KcMz/kB- This has the same form as the FDR in an 
equilibrium system except for T e ff which replaces the 
equilibrium temperature. 



IV. NUMERICAL SIMULATIONS 

In this section, we compare our three theoretical pre- 
dictions described in the previous section with results of 
numerical simulation of a two-dimensional granular gas 
system, which are the power spectrum Eq. ([8]). the fre- 
quency response function Eq. (|10p. and the fluctuation- 
dissipation relation Eq. (fT3| . The numerical setting of 
our system consists of N inelastic hard disks of mass m 
and diameter a moving in two-dimensional system on a 
thermal wall with a fixed temperature To. Here, x- and 
z-axes are set at horizontal and vertical directions of the 
system, respectively. The system width is denoted as L 
and periodic boundary conditions arc adopted in the hor- 
izontal direction at x = and x = L. The bottom wall is 
located at z = and there is no top wall. Gravitational 
force exerts on each disk along negative z-direction. In- 
clastic collisions between hard disks are considered by 
normal restitution coefficient r. When a disk collides 
with a thermal wall at the bottom, it comes off with a 
value of z-component of velocity v z sampled from the 
probability density 



P(v z ) 



fcsT 



exp 



2k B T Q 



(14) 



The horizontal component of velocity does not change 
during the collision with the thermal wall. 



Numerical simulations were performed with an event- 
driven algorithm devised to enhance speed of calcula- 
tion in particular in dense hard sphere systems [28|. In 
the following, all simulation data is presented with mass, 
length and time in units m, a, and a / [hsT® / m) 1 / 2 , re- 
spectively. This corresponds to choosing ksTo = 1. We 
set N = 5000 and L = 100, and keep these parame- 
ters unchanged throughout this paper. For our main 
results, r = 0.99 ~ 0.999 and g = 10~ 3 are used un- 
less otherwise mentioned. These correspond to N z = 50, 
0.05 < X < 0.5, and 1.98 < A < 6.25. A system with 
the width L = 100 and r > 0.99 is small enough to 
prevent any horizontal pattern structures (such as rip- 
ple and undulation). The global temperature T and the 
thermal velocity c are calculated using the following for- 
mula: T = K/ks and c = (2AT/m) 1 / 2 , where K is the 
long time average of the kinetic energy per disk. 



A. Macroscopic properties in NESS 

In Fig. [T] (top) , we show typical snapshots of particle 
configuration in the system of N = 5000 and g = 10~ 3 
for r = 0.999 and 0.992. The corresponding area-fraction 
profiles are plotted in Fig. Q] (bottom) . For a nearly elas- 
tic case (r = 0.999), the profile has one peak around the 
height z ~ 350, however, the area fraction is relatively di- 
lute less than 0.06 even at the height of the peak. Many 
inelastic particles are raised up relatively high like the 
equilibrium profile of the Boltzmann distribution. Con- 
trastingly, in case of r = 0.992, the profile drastically 
changes. The most of particles condense at relatively 
low height like a cluster; the area-fraction profile show a 
clear peak upon the low density region around the ther- 
mal wall. This state is known as density inversion state 
which has been observed in many experiments [29l . |30T | 
and simulations jU HH of vibrofluidizcd granular mat- 
ter. 

In accordance with the theoretical study by Bromberg 
ct al. UK showing that the steady state is characterized 
by a single parameter A defined in Eq. ([2"]l. we plotted 
in Fig. [5] the average kinetic energy per disk K and the 
average kinetic energy of the COM Kcmz as a function 
of A. The statistical error bars with the standard devia- 
tion were also plotted in all figures throughout the paper. 
A = (i.e., r = 1) corresponds to the equilibrium state 
where equipartition of energy IKcmz = K = ksTo = 1 
is satisfied. The factor 2 comes from the fact that Kqmz 
is defined using only z-component of COM velocity. Note 
that the horizontal component of velocity of the COM 
vanishes in our simulations because horizontal compo- 
nent of velocity of disks arc unchanged at collisions with 
the bottom wall. While K systematically decreases as a 
power law ~ A -148 , Kcmz shows a minimum at A ~ 2 
(r = 0.999) and increases with A for A > 2. Fig. [2] 
clearly indicates that equipartition of energy breaks down 
when A > 2 (i.e., IKcmz 7^ K). Similar behavior in 
much smaller systems has been preliminarily reported by 
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FIG. 1: Top: Snapshots of the two-dimensional simulation 
with N = 5000, L = 100 and g = 10~ 3 for different values of 
restitution coefficient (a) r = 0.999; (b) r = 0.992. Bottom: 
Area fraction profiles averaged over a long time period for 
r = 0.999 and r = 0.992. 



ours [2l[. In Ref. [23(, it has shown that the density in- 
version appears above the threshold A c (A > A c ) where 
A c ~ 1.06569. In the density inversion state, which be- 
comes pronounced for A > 2 as shown in Fig. Q] low 
density and high temperature gaseous region near the 
bottom may cause large fluctuations of the dense cluster 
on top of it. Therefore, this violation of the equipartition 
of energy should be closely connected to development of 
density inversion. 

Relation between the long time average of the COM 
height Z and the kinetic energy per particle is given in 
Fig. 3. A linear relation Z = K/mg + const, is well sat- 
isfied for K > 0.1 even if the system shows density inver- 
sion with a relatively high density cluster. In equilibrium 
system of dilute gas, the relation Z = K/mg + const, 
holds as a result of statistical mechanics. The fact that 
K characterizes Z in the same way as in equilibrium sug- 
gests that the global granular temperature T in inhomo- 
geneous nonequilibrium state still retains the same mean- 
ing as the equilibrium temperature at least in a macro- 
scopic sense. 
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FIG. 2: Kinetic energy per particle K (circles) and kinetic 
energy of the COM Kcmz (squares) are plotted versus A for 
r = 0.9999, 0.9996, 0.999, 0.998, 0.996, 0.994, 0.992, and 0.99 
from left to right. The solid line gives a numerical fit of the 
form 1.04 x A -1 ' 48 . 




FIG. 3: The average height of the center of mass Z versus K 
for r = 0.9999, 0.9996, 0.999, 0.998, 0.996, 0.994, 0.992, and 
0.99 from right to left. The error bars were smaller than the 
size of the marks. The solid line gives a linear fit with the 
slope (mg) -1 , where m = 1 and g — 10~ 3 . 



In Fig. @J we plot if as a function of the gravitational 
acceleration g. The dependence of K on g turned out to 
be rather week. We utilize this fact when we measure 
response function from simulations in an efficient way 
(see C in this section). 

In Fig.0 we plotted the probability distribution P(C) 
of the scaled COM velocity C = V z /(2K C Mz/M) 1/2 - 
The data are fitted sufficiently well by a Gaussian for 
all cases studied in this paper as is expected from the 
central limit theorem. This Gaussian property is consis- 
tent with our theory based on a linear Langevin equation 
with an additive Gaussian noise. 
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FIG. 4: Kinetic energy per particle K is plotted as a function 
of g. The error bars were smaller than the size of the marks. 




FIG. 5: Probability distribution of the scaled COM velocity 
C = V z /(2K C Mz/M) 1/2 . The solid line is Gaussian with 
unity dispersion. 



B. Power spectrum of the COM height 

We first test the theoretical prediction Eq. (jHJ on the 
power spectrum of the height of the COM. Using the 
relation Eq. (J6j, Eq. ((8j) can be rewritten as, 



S(lj) = S{Cbg/c)/ 



4 [ c\ 3 K C Mz 



9 J M 



f) 2 -W 2 +(jUO>)' 



(15) 



where Cj is scaled angular frequency defined by il> = ujcj g. 
This expression suggests that if we scale the power spec- 
trum and the angular frequency as in Eq. (|15[) . it shows 
a universal behavior independent of any system parame- 
ters. 
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In Fig. [6] (top), we plot the power spectrum S(lu) for 
different values of r. The two sharp peaks are observed; 
one is near zero angular frequency (lu — 0), the other 
one is at the angular frequency of macroscopic oscillation 
(lu = lu osc ) which increases as r decreases. The height of 
both these two peaks decreases as r decreases. We show 
in Fig. [S] (bottom) the scaled power spectrum S(Cu) that 
was obtained by scaling S(uj) in Fig. [5] (top) according 
to Eq. (|15[) using c and Kcmz calculated from simula- 
tion data. The theoretical prediction Eq. ([5]) with fitting 
numerical parameters jl — 0.50, = 1.7 is presented 
as a thick solid line. It is consistent with the results 
of simulations for a range of 0.99 < r < 0.994 in this 
region near the peak at Cj = Cj osc = Lu osc c/g, in which 
we expect our theory serves a description to hrst order 
approximation. Wc find large deviations from the the- 
oretical prediction in a region of Cj < Cj osc : the sharp 
peak near Cj = 0. As we will illustrate below this peak 
near Cj = can be associated with slow fluctuations of 
global granular temperature due to thermal conduction 
and collisional dissipation. Since T t herm/T OS c ~ N z ^> 1 
and Tdiss/Tosc ~ [N z (l — r 2 )]^ 1 > 1.0 for our simulations 
with N z = 50 and r > 0.99, the contribution of these two 
processes should appear at Cj < Cj osc . We also find for 
r > 0.998, the simulation data in Fig. [6] (bottom) devi- 
ate from our theory even in the region near the peak at 
Cj = Cj osc . These deviations near Cj osc might be attributed 
to the drastic change in density profiles shown in Fig. [1] as 
r is varied. Concerning our theory, the change in density 
profiles may affect the numerical coefficients and jl in 
Eq. ([B]). Furthermore, in the region lu < lu osc , the effect 
of fluctuations of global temperature mentioned above 
can be pronounced for r > 0.998 because both Ttherm 
and Tdiss become much larger than t osc , and hence the 
fluctuations have long life time. We admit that satisfac- 
tory explanation of these deviations for r > 0.998 has not 
been given yet. 

We show here a result of simulations suggesting that 
behavior of S(lu) in the region near to = can be well 
described by taking into account slow dynamics of K(t). 
Let us denote the slowly varying part of K(t) as K'(t) 
and suppose it fluctuates with a much longer time scale 
than t osc due to thermal conduction and collisional dis- 
sipation. Then, K'(t)/kB can be regarded as a time- 
dependent global granular temperature. Similarly, let 
Z'(t) denote the slowly varying part of Z(t) with the 
same time scale as K'(t). We assume here that in this 
long time scale, K'(t) and Z'(t) play the same role as 
their long time averages K and Z. That is, they sat- 
isfy the same linear relation as their long time averages 
observed in Fig. [3] Z'{t) = K'(t)/mg + const, with the 
same constant factor. If this is the case, the power spec- 
trum of SZ(t), S(lu), in the region near lu = should 
be given by the power spectrum of SK'(t)/mg, where 
8K'{t) = K'{t) — K. In Fig. [3 we show the power spec- 
trum of 5K(t)/mg, where 5K(t) = K(t) - K, and S(oj) 
for the case of r = 0.999 and 0.992. The figure shows 
that the curves around the peak in S(lu) near lu = and 



c«r=0.999 
Q-sr=0.998 




(0 



1.5i 
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r=0.999 
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r=0.998 






r=0.996 




A 


r=0.994 




<l 


r=0.992 






r=0.99 






FIG. 6: Top: The power spectrum for the COM height plot- 
ted versus angular frequency lu. Bottom: The scaled power 
spectrum for the COM height. The solid line depicts the 
theoretical prediction given in Eq. (|15p with ji — 0.50 and 
Q = 1.7. 



the peak in the power spectrum of SK(t)/mg near lu = 
are consistent. The consistency between the two curves 
is also observed for the other r values. This result indi- 
cates that the peak in S(lu) near lu = can be accounted 
for by slow dynamics of K (t) due to thermal conduction 
and collisional dissipation. It should be stressed that 
in our present theory fluctuations of granular tempera- 
ture in both space and time have been ignored and only 
the global granular temperature T is defined using the 
long time average of K(t). Further investigation is nec- 
essary to construct a theory that fully describes behavior 
of S(lu) taking into account the effect of slow fluctuations 
of granular temperature. 



C. Response functions 

Next, we test the theoretical prediction Eqs. (TTTT) and 
(|T2l on the frequency response functions for the COM. 
By scaling these functions in the same way as the power 
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FIG. 7: The power spectrum (PS) of SK(t) /mg and the power 
spectrum of 5Z(t), S(u), for (a) r = 0.999 and (b) r = 0.992. 



spectrum, we can derive universal equations, 



X'M = X '(ug/c)Mg 2 /c 2 = - 2 , (16) 



X» ee X "(G>g/c)Mg*/c* = ^ .(17) 

Measurement of the frequency response function by 
numerical simulations was performed using the following 
procedure. First, we prepared for a system in the station- 
ary state with a given N z , r, and g after sufficiently long 
time relaxation from the initial state of particles with 
randomly distributed positions and velocities. At t = 0, 
we exerted a small constant external force to all parti- 
cles in the direction of gravity and measured the height 
of the COM at t > 0; from this relaxation process of 
the COM, we can deduce a response function by a stan- 
dard procedure given in textbooks (see, e.g. Ref. [33|). 
That is, we measured a response function against a step 
functional external force. The frequency response func- 
tion is obtained as the Fourier transform of the response 
function. 

It is important to note that in the response of the COM 
height against a small but finite external force in our 
system, nonlinear effects resulting from change of time 
scales are not negligible. This can be seen from the fact 
that the relevant time scales shown in Eq. (JXJ) all depend 
on g, and that exerting a constant force in the direction 



of gravity is equivalent to changing g. Therefore, a lin- 
ear response can be defined only in the limit of small 
external force. This is a point that our Langevin-type 
theory is different from the well-known Langevin theory 
for a Brownian motion in a fixed harmonic potential, in 
which response of Brownian particle is linear against a 
finite external force. Consequently, we have to exert an 
external force that is much smaller than the gravitational 
force in our system in order to measure linear response 
of the COM height. Because fluctuation of the COM 
height of the 5000 particles is typically much larger than 
response against such a small constant force, we need to 
perform this measurement of the response function for a 
large number of systems with the same parameters N z , r, 
and g but different initial conditions and take an average 
of the response functions over all realizations. As shown 
later, in a case of a constant force that is 1% of gravi- 
tational force, we needed more than 10 realizations to 
obtain sufficient statistics for getting the clear response 
functions. This requires relatively long CPU time that 
impedes carrying out simulations with a wide range of 
parameters r, N z and g. 

We therefore optimized an efficient method by choosing 
appropriate parameter to approximately evaluate the re- 
sponse function from small number of realizations, which 
can be provided by the acceptable time in our computa- 
tional facilities. Suppose a system with gravitational field 
g is initially in a NESS and the gravitational acceleration 
is increased at t = from g to g + Ag. This is equiva- 
lent to exert a step functional external force —AIAg9(t) 
on the COM height, where 9(t) is the Heviside unit step 
function. Now let us define a function x(^3;5 + Ag) as 

x (t;g,g + Ag) = -^±/MAg, (18) 

where (•••)* represents the average taken over the en- 
semble of realization SZ at time t. This is a function of 
Ag in our system due to the nonlinear effects mentioned 
above; it would equal to the response function only if 
(SZ) t were linear in Ag. Let us denote the Fourier trans- 
form of x(t; 9, 9 + Aff) as x(w; g, g + Ag). According to 
Eq. ([5]), the frequency response function x( w >3) f° r the 
system in the stationary state with g is given by 

X(u;g) = lim v(w;s,s + Ag). (19) 

Ag->0 

We now consider the time scales that we introduced 
in Sec.IIB, which characterize macroscopic dynamics at 
t > 0. As we discussed in Sec.IIB, all these time scales in 
the NESS depend on g in the form r = c(g)/g x const., 
where we wrote g dependence of c explicitly for the sake 
of clarity. On the basis of our observation in Fig. 2] that 
c(g) changes a few percent as g is increased 10%, we 
assume thermal velocity at t > is given by c(g) if Ag is 
sufficiently small. Thus, these time scales at t > have 
the form r = c(g)/(g + Ag) x const., where g + Ag is 
gravitational acceleration at t > 0. This dependence of 
all the characteristic time scales on Ag leads us to the 
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following scaling relation: 



X(u;g,g + Ag) = 




where x is a nondimensional function. 

As far as the scaling relation Eq. (f^Uj) holds, we can 
estimate the limit in Eq. (jT9")l as follows: 

( g_+Ag \ 2 ( g + Ag \ 
X{u>9)=[ — - — ) X I w — - — ;g,g + Agj. (21) 

In order to verify validity of the scaling relation 
Eq. (|20]l. we carried out two series of simulations for r = 
0.992: Firstly, we measured the function x( w ! g-, 9 + Aff) 
in Eq. (|T5|) for Ag/g = 10~ 2 taking the average over 
41000 realizations. Secondly, we measured the frequency 
response function x(w; <?) using Eq. (f2"Tj) for Ag/g = 10 _1 
taking the average over 800 realizations. In Figs. we 
compare the frequency response functions obtained from 
these two series of simulations. We found that they are 
consistent although there are some discrepancies in x'( w ) 
near uj — 0. 

Another evidence of validity of the scaling relation is 
the fact that a FDR in an equilibrium system is well 
satisfied when we measured the frequency response func- 
tion using Eq. (|2"Tj) . which will be discussed later (see, 

Fig. nnj). 

The frequency response functions to be presented be- 
low were obtained using the scaling relation Eq. (|20]) (and 
Eq. I|2ip) by averaging over 800 realizations. In Figs. |H1 
we show the real (top) and imaginary (bottom) parts of 
the scaled response functions, x' and x", as functions 
of uj. Here values of c in Eq. (fTB")) and (JTTJ) are calcu- 
lated in a NESS without perturbation. The theoretical 
predictions Eqs. (fT6)) and (JTTJ) with the same (universal) 
fitting parameter as estimated in Fig. 6, p, = 0.50 and 
f2 = 1.7, are shown by thick lines. It appears that x(u>) 
is consistent with the theoretical predictions if r < 0.996. 

D. Fluctuation-dissipation relation 

To test FDR Eq. (fl~3]) predicted by our theory, we eval- 
uate the left and right hand sides independently using 
the results of simulations on S(uj) (Sec. II. B) and x"( w ) 
(Sec. II. C) presented in previous subsections. Note that 
Kcmz are measured in the NESS where S(uj) is mea- 
sured. 

At first, we confirmed that the FDR holds within the 
error bounds of the simulation result in the whole range 
of uj in the case of r — 1 in Fig. 1101 in which the stationary 
state is just the equilibrium state of elastic particles on 
a thermal wall. In Fig. [Til the left and right hand sides 
of the FDR are plotted as a function of uj for different 
r values. For all r (0.99 < r < 0.999), we found that 
the FDR holds within the error bounds in the range of 
higher frequency part of uj including a region near the 
peak uj ~ uj p , which is very close to uj osc defined as the 
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FIG. 8: The real part x' (top) and the imaginary part x" 
(bottom) of the frequency response function plotted versus 
angular frequency uj for r = 0.992. Circles are the data for 
Ag/g = 10~ 2 without using Eq. (|21[) ; averages have been ob- 
tained over 41000 realizations. Squares are for Ag/g = 10 _1 
using Eq. (|21|l with averages obtained over 800 realizations. 



angular frequency of a peak in S(oj). We stress here that 
we defined T e// as T eff = 2K C Mz/k B in the FDR. The 
quantitative agreement in Fig. [TTJ supports this definition 
of T e f f using Kcmz instead of using the global granular 
temperature T since T e ff is more than three times as 
large as T for r < 0.996 (see Fig.©. 

We found systematic deviations in a region uj < uj p . 
This deviations are attributed to an existence of a peak 
near uj = in S(uj) shown in Fig. [6j As we have discussed 
in Sec. IV B, the peak near uj — could be connected with 
slow fluctuations of granular temperature due to thermal 
conduction and collisional dissipation. For r = 0.999, 
time scales of these two processes, r t ^ erm and t<h S s be- 
come much larger than r osc , and thus the fluctuations 
with long life time might be responsible for the devia- 
tion at small uj. For r < 0.994, the deviation seems to 
increase as r decreases. Since the system has less gran- 
ular temperature for smaller r, larger heat current from 
the thermal wall is induced, and it could cause larger 
fluctuations in global granular temperature. Further in- 
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FIG. 9: The real (top) and imaginary part (bottom) of the 
scaled frequency response function plotted versus scaled an- 
gular frequency Cj = ujc/g. Averages have been taken over 
800 realizations. The thick lines are the theoretical predic- 
tion Eqs. (|16[1 and (|17[) with fitting parameters: fi = 0.50, 
tl = 1.7. 



vcstigation is necessary to understand more precisely this 
violation of the FDR. 



V. CONCLUSION 

In this paper we have studied validity of the 
fluctuation-dissipation relation with regard to the COM 
motion in a NESS of a driven granular fluid under gravity. 
By neglecting fluctuations of global temperature caused 
by thermal conduction and collisional dissipation, which 
change much slower than macroscopic oscillation of the 
fluid, we have derived a Langevin equation for the COM 
height. This Langevin equation predicts functional forms 
of the correlation and response functions for the COM 
height that contain two phenomenological numerical con- 
stants (X and Cl to be used as fitting parameters. It also 
gives a fluctuation-dissipation relation accompanied with 
an effective temperature T e ff that characterizes agitat- 
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FIG. 10: The left hand side uS(w) /2k B T ef f and the right 
hand side x"(w) of Eg. (fl3)l for r = 1. N = 5000, L = 100 
and g = 10~ 2 with averages over 400 realizations for S(uj) 
and over 800 realizations for x"( w )- 
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FIG. 11: The left hand side uS(u)/2k B T eff and the right 
hand side x"H of E q- G3} for ( a ) T = 0.999, (b) r = 0.998, 
(c) r = 0.996, (d) r = 0.994, (e) r = 0.992, and (f) r = 0.99. 
N = 5000, L = 100 and g = 10~ 3 with averages over 400 
realizations for S(lu) and over 800 realizations for x"(w). 
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ing motion of the COM height by T e ff = 2Kcmz/^b- 

In order to test the fluctuation-dissipation relation, we 
have performed event-driven MD simulations and mea- 
sured the power spectrum and response function for the 
COM height. While the power spectrum was consistent 
with our theory for r < 0.994 and lj around the angular 
frequency of the slowest oscillation of the COM, it also 
showed large deviations from the theoretical prediction 
near ui = and deviations for r > 0.994. The response 
function agreed well with our theory for r < 0.996 but 
showed deviations for r > 0.996. Furthermore, we com- 
pared the left and right hand sides of the fluctuation- 
dissipation relation. The results showed that they were 
consistent in a region of oj near the largest peak, but 
deviations occurred near uj = 0. These deviations are 
originated from the power spectrum for the COM height, 
showing a sharp peak near lo — 0, which cannot be de- 
scribed by our theory. 

We showed that these deviations near uj = can be 
attributed to slow fluctuations of global temperature de- 
fined as slowly varying part of kinetic energy per particle 
K(t) due to thermal conduction and collisional dissipa- 
tion; these fluctuations of global temperature have been 
neglected in our theory. The deviations in the power 
spectrum and resulting violation of the FDR are expected 
to be accounted for by a theory that describes both of 
Z(t) and Kit), which will be investigated in our future 
study. 

In Ref. [34| a formula that connects violation of the 
FDR in a NESS with the energy dissipation or equiva- 
lently energy input from outside has been proposed. A 
theory extended to include slow dynamics of K(t) and 
direct measurement of energy input in our simulations 
might give some insight to generality of their formula. 

Finally, we mention that a basic question whether the 
effective temperature T e ff obtained here has any phys- 
ical meaning as temperature in thermodynamics is still 
unclear. Definition of effective temperature in a system 
that relaxes in several time scales, typically glass, have 
been argued in Ref. @, HH]. It would be interesting to 
apply their theory to our problem with three time scales 
Ttherrm T~di S s and t osc . It would become also important to 
investigate via extensive simulation what happens if two 
systems with different effective temperature are in con- 
tact with each other. These studies would measure the 
direction of heat flow directly which might clarify an in- 
sight and physical meaning of the effective temperature. 
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Appendix A: Langevin equation 

In this section, we summarize derivation pol |2~H of the 
Langevin equation that describes the motion of the COM 
of grains. 

The equation of motion for the COM of grains in the 
model described in Sec. II can be written in the following 
form: 

d 2 Z 

M— = -AIg + F b . (Al) 

The right hand side of the equation of motion for the 
COM must be, in general, the sum of the external forces 
acting on the grains. In our model, they are the grav- 
itational force —Mg and the z-component of the force 
exerted by the bottom wall F b . Thus, it is essential to 
understand property of F b for the study of COM motion. 

Let us consider its reaction force F£(= —Fb), the force 
exerted by grains against the bottom wall. A snapshot of 
the granular fluid is sketched in Fig. [12] (a). Suppose the 
COM is at a height Z and is moving downward, for in- 
stance, with a velocity V. Let us now change the frame of 
reference to the center of mass frame (see Fig. [13(b)). It 
looks then the bottom wall that lies the distance Z away 
from the COM in the z-dircction, is moving upward with 
the velocity — V. Now the problem is how to determine 
the force acting on the bottom wall as a result of 
frequent collisions of granular particles, in the situation 
shown in Fig. [12] (b), in which the bottom wall is moving 
upward with the velocity — V against the macroscopically 
static fluid. 




FIG. 12: (a) a sketch of the system observed in the labora- 
tory frame of reference, (b) the same system observed in the 
center of mass frame, (c) the Rayleigh piston: a piston un- 
dergoes random collisions with a one-dimensional heat bath 
of particles. 

This problem is similar to the problem of determining 
the force acting on a one-dimensional Brownian particle, 
which is referred to as the Rayleigh piston (l9| . moving 
in the z-directions with a velocity — V (sec Fig. [T2] (c)). 
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We model an expression of F' h on the basis of this analogy 
with Brownian motion, and assume that F£ consists of 
the following three components: The first is a systematic 
force fp(t) which equals to the pressure multiplied by 
the area of the bottom wall. Since local density near 
the bottom wall changes according to the motion of the 
COM, this force may depend on time. Apparently the 
long time average of fp(t), say fp, must be equal to 
Mg, the gravitational force acting on all particles. The 
simplest assumption on the time dependent part of fp(t) 
is that it is proportional to deviation of the COM height 
from its stationary value, Z(t) — Z, since change in the 
local density near the bottom wall could be proportional 
to -(Z(t) - ~Z) if the change Jn the height of the COM 
is sufficiently small: \Z(t) — Z\/Z <C 1. The second is 
a frictional force. We assume here the simplest form of 
the frictional force that is linear in the relative velocity 
-V{t) of the bottom wall to the COM. The third is a 
random force. Thus, we assume the following form: 



Fl(t) = -Mg + Mn 2 
= -F b (t), 



Z(t)- Z) +MfxV{t)+R'(t) 

(A2) 



where il is a coefficient that specifies the angular fre- 
quency of the slowest oscillation of the COM, and /i is 
the frictional coefficient. According to the discussion on 
the characteristic time scales in Sec.IIB, the time scales 
for macroscopic oscillation t osc and that for pressure re- 
laxation t p are t osc ~ r p ~ c/g. Thus, we assume 

fl = &/t osc = Clg/c, /i = £/t p = fxg/c. (A3) 

As property of the random force, we assume stationary 
Gaussian white noise in the same way as the Rayleigh 
piston: 

(R'(t)) = 0, (R'{t)R'(t')) = IS(t - t'). (A4) 

where / represents the intensity of the random force. 

Substituting F^ obtained in (|A2[) into the equation of 
motion of the COM ([ATI) , we obtain 



d 2 SZ 
dt 2 



-Q 2 6Z 



/'- 



dbZ R(t) 



dt 



M 



(A5) 



where SZ = Z(t) - Z and R(t) = -R'(t). The random 
force R{t) has the exactly same property described in 
(|A4[) as R'(t). Note that the Langevin equation (|A5[) has 
the same form as that describes Brownian motion in a 
harmonic potential. 



Appendix B: Derivation of the power spectrum and 
the response function 

Since derivation of the power spectrum and the re- 
sponse function from the Langevin equation describing 
Brownian motion in a harmonic potential is given in text- 
books (see, e.g., Ref. [36|), we present here only essential 
steps in their calculation. Let us first consider the power 
spectrum of the fluctuating motion of the COM obey- 
ing the Langevin equation (|A5|) . The formal solution of 



Eq. (|A5|) is written as 

Z(t)-Z= j G(t-t')^p-dt' + F ini (t), (Bl) 



where the function G(t) is given by 



G(t.) = 



w 



sin (iv t) , 



(B2) 



and u is defined by uj q = (ft 2 - (fi/2) 2 ) 1 / 2 . The last 
term F ini (t) in Eq. (|B1|) consists of those that depend on 
the initial conditions and vanish after a sufficient amount 
of time. Thus, the term is negligible when calculating 
long-time averages of physical quantities in the stationary 
state. 

Using this formal solution, one can calculate the two- 
time correlation function <f>(t) in a NESS defined by 
<j)(t) = lim t '-).oo (6Z(t')SZ(t' + t)), where the brackets 
(• ■ ■ ) indicate an average over the random force R(t). 
We took the limit t' — > oo in order to assure the system 
to be in the stationary state. 

The power spectrum of SZ(t) can be obtained using 
the Winner-Khinchin theorem: 



S(cj) 



dte- iut (j){t) 
I 1 



M 2 (02 _ ^2)2 + (ajw) 2 



(B3) 
(B4) 



Next, let us consider the response function for the 
COM, which describes linear response of the COM with 
regard to a small external force sf(t). The Langevin 
equation in this case is written as 

d 2 6Z 2 d5Z sf(t) R(t) 

Taking the average over the random force, we obtain 



The response function x(t) is defined as 



(SZ(t)} 



dt' x (t - t')ef(t'). 



= 0. (B6) 



(B7) 



Here the external force ef(t) is assumed to be infinitely 
small. The Fourier transform of this relation yields 
(8Z(uj)) = x(w)e/(u;), and hence 



(B8) 



X (w) = M{iZ(«))/£/(w). 



Performing the Fourier transform of the relation (|B6[) 
and comparing it with Eq. (|B8|) , we obtain the frequency 
response function (complex admittance) 



1 



1 



M ft 2 



IfJIO 



(B9) 
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